/* Simple C implementation of the mandlebrot set */

#include <stdio.h>
#include <stdlib.h>
#include "print_png.h"
#include <string.h>

// Maximum number of iterations
#define MAX_ITER 255

extern double wtime();

int main (int argc, char* argv[]){
  
  // Will be used to create the image
  char *values; 

  // Number of pixels in the generated square
  int size;
 
  // default values
  long double xMin = -1.5;
  long double xMax = 0.5;
  long double yMin = -1.0;
  long double yMax = 1.0;

  double t1, t2, t3, diff, print_time;

  char fileName[50];

  if((argc!=3)&&(argc!=7)){
    printf("Usage: mandbrot.x $SIZE $FILENAME [$XMIN $XMAX $YMIN $YMAX]\n");
    printf("If mins and maxes not entered they will default to\n");
    printf("-1.5, 0.5, -1.0, and 1.0 respectively\n");
    exit(1);
  }
  
  size=atoi(argv[1]);
  
  strcpy(fileName,argv[2]);

  if(argc==7){
    xMin=atof(argv[3]);
    xMax=atof(argv[4]);
    yMin=atof(argv[5]);
    yMax=atof(argv[6]);
  }
  
  values = (char*) malloc (sizeof(char)*size*size);
  if(values==NULL){
    fprintf(stderr,"Malloc failed\n");
    exit(1);
  }

  t1=wtime();
  
  calc_mb(xMin, xMax, yMin, yMax, values, size);

  t2=wtime();

  print_png(values, size, size, fileName);

  t3=wtime();

  diff=t2-t1;
  print_time=t3-t2;

  printf("Time to run program = %f\n",diff);
  printf("Time to print image = %f\n",print_time);
  printf("%dx%d matrix\n",size,size);
  printf("xmin = %Lf, xmax = %Lf\n",xMin,xMax);
  printf("ymin = %Lf, ymax = %Lf\n",yMin,yMax);

  return 0;
}

/* Function calc_mb
   Calculates the number of iterations to determine whether or not a point
   is a part of the mandlebrot set
   Args xMin, xMax, yMin, yMax - map out the area for which the mandlebrot
     set is being calculated.  Not changed.
   Arg values - array in which the calculated values will be stored.
     Values are overwritten on output.
   Arg size - size of the nxn values array.  Not changed.
 */
int calc_mb(long double xMin, long double xMax, long double yMin,
	    long double yMax, char *values, int size){

  int i, j, k;
  long double x, y;
  long double rn0, rn1, sn0, sn1;
 
 
  for(i=0; i<size; i++){
    for(j=0; j<size; j++){
      x = ((long double)i/size)*(xMax-xMin) + xMin;
      y = ((long double)j/size)*(yMax-yMin) + yMin;
      rn0 = x;
      sn0 = y;
      for(k=0; k<MAX_ITER; k++){
	
	rn1=rn0*rn0-sn0*sn0+x;
	sn1=2*rn0*sn0+y;

	//printf("x = %6.3Lf, y = %6.3Lf, rn1 = %6.3Lf, sn1 = %6.3Lf\n",x,y,rn1,sn1);

	// rule of thumb cutoff that shows a point will not be in the set
	if((rn1*rn1+sn1*sn1)>4.0){
	  break;
	}
	
	rn0=rn1;
	sn0=sn1;


	
      }
      values[j*size+i]=k;
    }
  }
}
